Problem statement

Overview

explosion.png

A gas canister explodes on a flat surface. Several small pieces of debris are sent flying in all directions. Where will they land?

Modelling assumptions

Useful formulas

Finding the distribution for the range $ x $

From the definition of the CDF, we have $ F_X(x) = P(X \leq x) = P(g(U, \Theta) \leq x) $. To work out the RHS, we integrate the joint PDF over the region of all $ u $ and $ \theta $ which gives $ X $ less than $ x $:

$$ F_X(x) = \int \limits_{(u, \ \theta): \ g(u, \ \theta) \ \leq \ x} f_{u, \ \theta}(u, \ \theta) \ \text{d} u \ \text{d} \theta. $$

where $ f_{u, \ \theta}(u, \ \theta) $ is the joint probability density function, which depends on the distributions we choose to model the random variables $ u $ and $ \theta $ as having. We will leave this as an arbitrary function for now.

In order to work out the region of integration, we will look at the contour plot for $ g(u, \ \theta) $.

After seeing the contour plot, it will be easier to integrate over the region where $ x $ is larger than some value, then subtract it from 1. So our formula will use

$$ F_X(x) = 1 - \int \limits_{(u, \theta): \ g(u, \theta) \geq x} f_{u, \theta}(u, \theta) \ \text{d} u \ \text{d} \theta $$

We now need to find the limits of integration. On the $ u $-axis, the region will be from the left-most peak of the contour, all the way up to infinity. This peak will always have $ \theta = \frac{pi}{4} $, which allows us to easily solve for value of $ u $ at the peak:

$$ x = \frac{u^2 \sin \left (2 \times \frac{\pi}{4} \right )}{g} \ \ \rightarrow \ \ u = \sqrt{gx}. $$

Next, for a given value of $ u $, we need the range of $ \theta $. This will be a closed interval between the lower and upper intersections with the contour. The lower limit and upper limit will be, respectively,

$$ \theta_{-} = \frac{1}{2} \sin^{-1} \frac{gx}{u^2} \ \ \text{and} \ \ \theta_{+} = \frac{1}{2} \cos^{-1} \frac{gx}{u^2}. $$

So our formula is

$$ F_X(x) = 1 - \int \limits_{\sqrt{gx}}^{\infty} \int \limits_{\frac{1}{2} \sin^{-1} \frac{gx}{u^2}}^{\frac{\pi}{2} - \frac{1}{2} \sin^{-1} \frac{gx}{u^2}} f_{u, \ \theta}(u, \ \theta) \ \text{d} \theta \ \text{d} u. $$

If we know $ f_{u, \ \theta}(u, \ \theta) $, then this is a formula which we can work with...

Choosing our Statistical Model

There is lots of room for experimentation here i.e. how does changing the input distributions affect the distribution of projectile landing points. In this particular model, we will use:

Then we have the following marginal PDFs:

$$ f_U(u) = \begin{cases} \frac{u}{144} \exp \left ( -\frac{u^2}{288} \right ), & u \geq 0, \\ 0, & \text{otherwise} \end{cases} \ \ \text{and} \ \ f_{\Theta}(\theta) = \begin{cases} \frac{4 \pi - |4 \pi - 16 \theta |}{\pi ^2}, & 0 \leq \theta \leq \pi/2, \\ 0, & \text{otherwise} \end{cases} $$

Since $ u $ and $ \theta $ are independent, the joint PDF is simply their product:

$$ f_{U, \ \Theta} = f_U(u) \times f_{\Theta}(\theta) = \begin{cases} \frac{u}{144} \exp \left ( -\frac{u^2}{288} \right ) \frac{4 \pi - |4 \pi - 16 \theta |}{\pi ^2}, & u \geq 0 \ \ \text{and} \ \ 0 \leq \theta \leq \pi/2, \\ 0, & \text{otherwise} \end{cases}. $$

These functions are implemented and plotted here.

We can now implement the CDF of $ x $:

Plot the CDF of $ x $:

Next, the PDF. Instead of differentiating this numerically, which would be slow and inaccurate due to the complexity of the function, we will manipulate algebraically first. We have

$$ f_X(x) = \frac{\mathrm{d} }{\mathrm{d} x} F_X(x) = - \frac{\mathrm{d} }{\mathrm{d} x} \int \limits_{\sqrt{gx}}^{\infty} \int \limits_{\frac{1}{2} \sin^{-1} \frac{gx}{u^2}}^{\frac{\pi}{2} - \frac{1}{2} \sin^{-1} \frac{gx}{u^2}} f_{u, \ \theta}(u, \ \theta) \ \text{d} \theta \ \text{d} u. $$

We will need to use the Leibniz rule twice to expand this. First, treat the inner integral as an arbitrary function $ F(x, u) $:

$$ \frac{\mathrm{d} }{\mathrm{d} x} \int \limits_{\sqrt{gx}}^{\infty} F(x, u) \ \text{d} u = -F(x, \sqrt{gx}) \times \frac{g}{2\sqrt{gx}} + \int \limits_{\sqrt{gx}}^{\infty} \frac{\partial F}{\partial x} \ \text{d} u $$

Now, work out the partial derivative of $ F $:

$$ \frac{\partial F}{\partial x} = \frac{\partial }{\partial x} \int \limits_{\frac{1}{2} \sin^{-1} \frac{gx}{u^2}}^{\frac{\pi}{2} - \frac{1}{2} \sin^{-1} \frac{gx}{u^2}} f_{u, \ \theta}(u, \ \theta) \ \text{d} \theta = -\frac{g}{2\sqrt{u^4 - g^2 x^2}} \left ( f_{u, \ \theta}(u, \ \frac{\pi}{2} - \frac{1}{2} \sin^{-1} \frac{gx}{u^2}) + f_{u, \ \theta}(u, \ \frac{1}{2} \sin^{-1} \frac{gx}{u^2}) \right ) $$

Since $ F(x, \sqrt{gx}) $ = 0, we get our PDF as:

$$ f_X(x) = \int \limits_{\sqrt{gx}}^{\infty} \frac{g}{2\sqrt{u^4 - g^2 x^2}} \left ( f_{u, \ \theta}(u, \ \frac{\pi}{2} - \frac{1}{2} \sin^{-1} \frac{gx}{u^2}) + f_{u, \ \theta}(u, \ \frac{1}{2} \sin^{-1} \frac{gx}{u^2})\right ) \ \text{d} u. $$

This function is implemented now.

Now we can plot the PDF of $ x $:

At first glance this seems reasonable, but let's evaluate some summary statistics.

Finding some properties of the distribution of $ x $

Let's check more carefully, by looking at the mean range, $ E[X] $. By LOTUS, we should get

$$ E[X] = \int \limits_{-\infty}^{\infty} \int \limits_{-\infty}^{\infty} g(u, \theta) \times f_{u, \ \theta}(u, \ \theta) \ \text{d} u \ \text{d} \theta $$

while from our PDF, we should get

$$ E[X] = \int \limits_{-\infty}^{\infty} x f_X(x) \ \text{d} x. $$

We can calculate both of these values explicitly and check they are equal:

We can see that they are in fact equal - to within 3 d.p.

We can also get the variance:

By visual inspection of the graph of the PDF again, these values seem reasonable.

Next we can look at the median and quartiles. This will be easiest done by reading from the CDF.

Since the median is less than the mean (15.58 < 23.80), the distribution is skewed right (positive skew), which again matches what we see in the PDF graph.

Finally, calculate skewness and kurtosis quantitatively.

This confirms the positive skew, and we also see a positive excess kurtosis (leptokurtic), which means the tails are less filled out relative to the Normal distribution.

Simulation with Random Numbers (Monte-Carlo)

We can also simulate the problem directly, by sampling from the chosen distributions and calculating range for a very large number of trials. We will use 100,000 trials here, and also add a degree of freedom to give a radially-symmetric distribution. The top-down angle (written as $ \phi $) will be uniformly distributed on $ [0, 2 \pi) $.

Finally, let's plot our simulated landing distribution:

The densely-populated dark blue ring close to the centre shows the region where most of the pieces land. This is consistent with our distribution being more populated at smaller ranges. The marginal distributions, which should match our PDF, also resemble the correct shape.